###差异分析高低风险组，火山图


#引用包
library(limma)
expFile="OV_TPM.txt"          #表达数据文件
riskFile="rs_tcga.txt"       #风险文件
logFCfilter=1                 #logFC过滤条件
fdrFilter=0.05                #fdr过滤条件
setwd("E:\\1.Ovarian_Treg\\1.picture\\16.go")    #设置工作目录

#读取表达数据文件,并对输入文件整理
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp), colnames(exp))
data=matrix(as.numeric(as.matrix(exp)), nrow=nrow(exp), dimnames=dimnames)
data=avereps(data)

##去除正常样品
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
data=t(data)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(data))
data=avereps(data)
data=t(data)

#读取risk文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)
risk$risk=ifelse(risk$riskScore>median(risk$riskScore),"high","low")
sameSample=intersect(colnames(data), row.names(risk))
data=data[,sameSample]
risk=risk[sameSample,]
#colnames(risk)[13]="risk"
#提取低风险组和高风险组的样品
riskLow=risk[risk$risk=="low",]
riskHigh=risk[risk$risk=="high",]
dataLow=data[,row.names(riskLow)]
dataHigh=data[,row.names(riskHigh)]
data=cbind(dataLow,dataHigh)
data=data[rowMeans(data)>1,]
conNum=ncol(dataLow)
treatNum=ncol(dataHigh)
Type=c(rep(1,conNum), rep(2,treatNum))

#差异分析
outTab=data.frame()
for(i in row.names(data)){
  rt=data.frame(expression=data[i,], Type=Type)
  wilcoxTest=wilcox.test(expression ~ Type, data=rt)
  pvalue=wilcoxTest$p.value
  conGeneMeans=mean(data[i,1:conNum])
  treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])
  logFC=log2(treatGeneMeans)-log2(conGeneMeans)
  conMed=median(data[i,1:conNum])
  treatMed=median(data[i,(conNum+1):ncol(data)])
  diffMed=treatMed-conMed
  if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){  
    outTab=rbind(outTab,cbind(gene=i,lowMean=conGeneMeans,highMean=treatGeneMeans,logFC=logFC,pValue=pvalue))
  }
}
pValue=outTab[,"pValue"]
fdr=p.adjust(as.numeric(as.vector(pValue)), method="fdr")
outTab=cbind(outTab, fdr=fdr)
save(outTab,file = "all.diff.rdata")
write.table(outTab,file = "all.diff.txt",row.names = F,quote = F,col.names = T,sep="\t")
#输出差异表格
outDiff=outTab[( abs(as.numeric(as.vector(outTab$logFC)))>logFCfilter & as.numeric(as.vector(outTab$fdr))<fdrFilter),]
write.table(outDiff, file="riskDiff.txt", sep="\t", row.names=F, quote=F)


#渐变火山图
# 加载包：
library(ggplot2)
library(ggrepel)
# 读取数据：
data=read.table('all.diff.txt',sep = "\t",header = T)
model_gene=read.table('modelgene.txt',sep = "\t",header = T)[,1]
# 加载包：
library(ggplot2)
index=order(data$pValue)
data=data[index,]
data$label <- c(model_gene,rep(NA,(nrow(data)-4)))

colnames(data)
ggplot(data,aes(logFC, -log10(fdr)))+
  # 横向水平参考线：
  geom_hline(yintercept = -log10(0.0001), linetype = "dashed", color = "#999999")+
  # 纵向垂直参考线：
  geom_vline(xintercept = c(-1.2,1.2), linetype = "dashed", color = "#999999")+
  # 散点图:
  geom_point(aes(size=-log10(fdr), color= -log10(fdr)))+
  # 指定颜色渐变模式：
  scale_color_gradientn(values = seq(0,1,0.2),
                        colors = c("#39489f","#39bbec","#f9ed36","#f38466","#b81f25"))+
  # 指定散点大小渐变模式：
  scale_size_continuous(range = c(0.5,3))+
  # 主题调整：
  theme_bw()+
  # 调整主题和图例位置：
  theme(panel.grid = element_blank(),
        legend.position = c(0.1,0.8),
        legend.justification = c(0,1),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7)
  )+
  # 设置部分图例不显示：
  guides(col = guide_colourbar(title = "-Log10_q-value"),
         size = "none")+
  geom_text_repel(aes(label = label),
                  size=2,color="#3C5488FF",
                  direction="both",
                  min.segment.length = 0.05,face="bold.italic",
                  segment.alpha=0.6,
                  label.padding = 0.4,
                  max.overlaps =50,nudge_x = 0.2,
                  nudge_y=0.2)+
  # 添加标签：
  #geom_text_repel(aes(label=label, color = -log10(fdr)), size = 3, vjust = 1.5, hjust=1)+
  # 修改坐标轴：
  xlab("Log2FC")+
  ylab("-Log10(FDR q-value)")

# 保存图片：
ggsave("vocanol_plot.pdf", height = 4, width = 4)
###go富集分析
##go富集分析


#引用包
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)
library(circlize)
library(RColorBrewer)
library(dplyr)
library("ggpubr")
library(ComplexHeatmap)

pvalueFilter=0.05       #p值过滤条件
qvalueFilter=0.05       #矫正后的p值过滤条件

#定义颜色
colorSel="qvalue"
if(qvalueFilter>0.05){
  colorSel="pvalue"
}

rt=read.table("riskDiff.txt", header=T, sep="\t", check.names=F)     #读取输入文件

#基因名字转换为基因id
genes=unique(as.vector(rt[,1]))
entrezIDs=mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
entrezIDs=as.character(entrezIDs)
gene=entrezIDs[entrezIDs!="NA"]        #去除基因id为NA的基因
#gene=gsub("c\\(\"(\\d+)\".*", "\\1", gene)

#GO富集分析
kk=enrichGO(gene=gene, OrgDb=org.Hs.eg.db, pvalueCutoff=1, qvalueCutoff=1, ont="all", readable=T)
GO=as.data.frame(kk)
#保存富集结果
GOBP=GO[GO$ONTOLOGY=="BP",]
GOBP=GOBP[order(GOBP$Count,decreasing = T),]
GOBP=GOBP[1:5,]
GOBP$Description[3]="	adaptive immune response"
GOCC=GO[GO$ONTOLOGY=="CC",]
GOCC=GOCC[order(GOCC$Count,decreasing = T),]
GOCC=GOCC[1:3,]
GOMF=GO[GO$ONTOLOGY=="MF",]
GOMF=GOMF[order(GOMF$Count,decreasing = T),]
GOMF=GOMF[1:3,]
GO=rbind(GOBP,GOCC)
GO=rbind(GO,GOMF)
colnames(GO)
GO$ONTOLOGY=factor(GO$ONTOLOGY,levels = c("BP","CC","MF"))
#载入ggplot2包；
library(ggplot2)
#install.packages("ggpubr")
#载入ggpubr包；
library(ggpubr)
#自定义图表主题，对图表做精细调整；
top.mar=0.2
right.mar=0.2
bottom.mar=0.2
left.mar=0.2

#隐藏坐标轴轴，并对字体样式、颜色、刻度长度等进行限定；
mytheme1<-theme_classic()+
  theme(text=element_text(family = "sans",colour ="gray30",size = 12),
        axis.line = element_blank(),
        axis.ticks = element_line(size = 0.6,colour = "gray30"),
        axis.ticks.length = unit(1.5,units = "mm"),
        plot.margin=unit(x=c(top.mar,right.mar,bottom.mar,left.mar),
                         units="inches"))

mycolor <- c("#0077c1","#00a99e","#6bc72b","#ff5a20","#ff1620","#752995")
p4 <- ggbarplot(GO, x = "Description", y = "Count",
                fill = "ONTOLOGY",
                color = "white",
                width = 0.9,
                palette = mycolor,
                sort.val = "asc",
                sort.by.groups = TRUE,
                label = TRUE,
                lab.size = 2.5,
                lab.col = "white",
                lab.vjust=0.5,lab.hjust=1.5,
                xlab = "",
                ylab = "Number",
                x.text.angle = 90,
                ggtheme = mytheme1)+coord_flip()
p4
ggsave("barplot.pdf",width = 8,height = 5)
